A partir d'une présentation de @AustinRochford¶
Inférence Bayésienne et Analyses fréquentistes¶
|
|
| Bayésienne | Fréquentistes |
Problème avec les p-values¶
- Gelman, A., & Geurts, H. M. (2017). The statistical crisis in science: How is it relevant to clinical neuropsychology?. The Clinical Neuropsychologist, 31(6-7), 1000-1014.
- Baker, M. (2016). 1,500 scientists lift the lid on reproducibility. Nature News, 533(7604), 452.
- Gelman, A. (2016). The problems with p-values are not just with p-values. The American Statistician, 70.
- Button, K. S., Ioannidis, J. P., Mokrysz, C., Nosek, B. A., Flint, J., Robinson, E. S., & Munafò, M. R. (2013). Power failure: why small sample size undermines the reliability of neuroscience. Nature - Reviews Neuroscience, 14(5), 365.
- Simmons, J., Nelson, L., & Simonsohn, U. (2011). False-positive psychology: Undisclosed fexibility in data collection and analysis allow presenting anything as signifcant. Psychological Science, 22, 1359–1366.
- Schatz, P., Jay, K. A., McComb, J., & McLaughlin, J. R. (2005). Misuse of statistical tests in archives of clinical neuropsychology publications. Archives of Clinical Neuropsychology, 20, 1053–1059
- Gelman, A. (2005). Analysis of variance—why it is more important than ever. The annals of statistics, 33(1), 1-53.
- Bezeau, S., & Graves, R. (2001). Statistical power and effect sizes of clinical neuropsychology research. Journal of Clinical and Experimental Neuropsychology, 23(3), 399-406.
- Meehl, P. E. (1990). Why summaries of research on psychological theories are often uninterpretable. Psychological Reports, 66, 195–244
Inférence Bayesienne¶
Exemple classique:
Supposons un test sanguin permettant de détecter le vampirimse avec un taux de succès de 95% Pr(positive|vampire) = 0.95. C'est donc un test précis, mais qui implique son lot de faux-positif. Un pourcent du temps le test sera positif pour une personne normale, donc Pr(positive|mortel)=0.01. La dernière donnée connue est qu'une personne sur mille est un vampire dans la population générale, donc Pr(vampire)=0.001
Probabilité d'être un vampire si on test positif? 95%? 1%?
Probabilité conditionnelle¶
En théorie des probabilités, une probabilité conditionnelle est la probabilité d'un événement sachant qu'un autre événement a eu lieu.
$$ \begin{align*} P(A\ |\ B) & = \textrm{La probabilité que } A \textrm{ arrive si on sait que } B \textrm{ est arrivé} \\ & = \frac{P(A \textrm{ et } B)}{P(B)}. \end{align*} $$
Notre question,
Quelle est la probabilité qu'une personne soit un vampire si elle est testée positive?
$$ \begin{align*} P(vampire) & = 0.001 \\ P(\textrm{Test } +\ |\ vampire) & = 0.95 \\ P(\textrm{Test } +\ |\ mortel) & = 0.01 \\ \\ P(vampire\ |\ \textrm{Test } +) & =\ \textbf{?} \end{align*} $$
Théorem de Bayes¶
Le théorem de Bayes nous montre comment passer de $P(B\ |\ A)$ à $P(A\ |\ B)$.
$$\large P(A\ |\ B) = \frac{P(B\ |\ A)\ P(A)}{P(B)}$$
Le terme P(A) est la probabilité a priori de A. Elle est « antérieure » au sens qu’elle précède toute information sur B. Le terme P(A|B) est appelée la probabilité conditionnelle de A sachant B . Elle est « postérieure », au sens qu’elle dépend directement de B. Le terme P(B|A), pour un B connu, est appelé la fonction de vraisemblance de A. De même, le terme P(B) est appelé la probabilité marginale ou a priori de B.
$$\color{red}{P(A\ |\ B)} = \frac{\color{blue}{P(B\ |\ A)}\ \color{orange}{P(A)}}{\color{green}{P(B)}}.$$
Pour plusieurs modèles la la probabilité marginale est impossible ou difficilement calculable analytiquement.
Pour le cas qui nous intéresse, on peut calculer analytiquement la probabilité avec le théorème de Bayes et la Formule des probabilités totales
$$Pr(vampire|+)=\frac{Pr(+|vampire)Pr(vampire)}{Pr(+)}$$
$$ Pr(+)=Pr(+|vampire)Pr(vampire)+Pr(+|mortel)(1−Pr(vampire)) $$
$$ \begin{align*} P(vampire) & = 0.001 \\ P(+\ |\ vampire) & = 0.95 \\ P(+\ |\ mortel) & = 0.01 \\ \\ P(vampire\ |\ \textrm{Test } +) & = \frac{0.95*0.001}{0.95*0.001 + 0.01*(1-0.001)} \end{align*} $$
(0.95*0.001)/(0.95*0.001 + 0.01 * (1-0.001))
C'est un exemple classique, mais en soi pas particuilièrement Bayésien¶
Il existe une façon plus intuitive de poser problème
Dans une population de 100,000 personnes, 100 sont vampires.
De ces 100 vampires, 95 vont tester positif au test.
Des 99,900 mortels, 999 vont tester positif au test.
Donc si on teste les 100,000 personnes, quelle est la proportion de ceux qui teste positif d'être des vampires? La question posée ainsi semble plus facile à répondre.
On a juste à compter le nombre de test positifs: 95 + 999 = 1094. On sait que dans le lot 95 sont des vampires:
Pr(vampire|+)=95/1094≈0.087
Programmation probabilistique pour inférence bayesienne¶
L'exemple précédant était déjà complexe à exprimer et pouvait se résoudre analytiquement parce qu'il était très simple.
Résolvons le même problème avec Python¶
On connait la prévalence de la maladie dans la population générale, soit 1/1000
La distribution de Bernoulli donne la probabilité binaire (pile ou face biaisé). Si $X \sim \textrm{Bernoulli}(p),$
$$ \begin{align*} P(X = 1) & = p \\ P(X = 0) & = 1 - p. \end{align*} $$
import pymc3 as pm
with pm.Model() as disease_model:
has_disease = pm.Bernoulli('est_vampire', 0.001)
Si la personne est vampire, il y a 95% de chance de tester positif. Sinon, il y a quand même 1% de chance de tester positif
with disease_model:
p_test_pos = has_disease * 0.95 + (1 - has_disease) * 0.01
La personne dans l'exemple a testé positive alors:
with disease_model:
test_pos = pm.Bernoulli('test_pos', p_test_pos, observed=1)
Quelle est la probabilité que cette personne ait la maladie
with disease_model:
disease_trace = pm.sample(draws=10000, random_seed=SEED, cores=1)
disease_trace['est_vampire']
disease_trace['est_vampire'].mean()
Ce qu'on a calculé tantôt:
(0.95*0.001)/(0.95*0.001 + 0.01 * (1-0.001))
Méthodes de Monte Carlo¶
Les méthodes de Monte Carlo utilisent un échantillionnage aléatoire pour aproximer des quantités difficile a obtenir analytiquement. Ils sont un outil essentiel pour l'inférence Bayesienne où la distribution marginale est habituellement difficile à obtenir.
Un exemple¶
On génère 5000 points distribués aléatoirement dans un carré unitaire
N = 5000
x, y = np.random.uniform(0, 1, size=(2, N))
fig
En comptant le nombre de point qui tombent dans le quatier de cercle de rayon 1, on obtient une approximation de l'aire dans de quartier de cercle, qui est en réalité de $\frac{\pi}{4}$.
in_circle = x**2 + y**2 <= 1
fig
4 * in_circle.mean()
Historique des méthodes de Monte Carlo¶
Les méthodes de Monte Carlo ont été largement utilisées pour le Projet Manhattan. Sur la photo ci-dessus, on voit des scientifiques du projet Manhattan Stanislaw_Ulam,Richard Feynman et John von Neumann. En travaillant sur le projet Manhattan, Ulam a donné l'une des premières descriptions des algorithmes des Markov Chain Monte Carlo. La librairie de programmation probabiliste Stan est nommé en son honneur.
Étude de cas : Le problème de Monty Hall¶
Le problème de Monty Hall est un célèbre puzzle de probabilité, basé sur l'émission de jeux des années 1960 Let's Make a Deal et nommé d'après son hôte d'origine. Dans le jeu, un concurrent se voyait présenter trois portes, dont deux contenaient un article de peu ou pas de valeur (par exemple, une chèvre) et l'une d'entre elles contenait un article de très grande valeur (par exemple, une voiture de luxe). Le concurrent devait d'abord deviner quelle porte contenait la voiture de sport. Après la devinette initiale du concurrent, Monty ouvrait l'une des deux autres portes, révélant une chèvre. Monty offrait alors au concurrent la chance de changer leur choix de porte. Le problème de Monty Hall pose la question suivante : le candidat doit-il conserver son choix initial de porte ou le changer ?
Au départ, nous n'avons aucune information sur la porte derrière laquelle se trouve le prix.
with pm.Model() as monty_model:
prize = pm.DiscreteUniform('prize', 0, 2)
Si nous choisisson la porte 1:
| Prix derrière | Porte 1 | Porte 2 | Porte 3 |
|---|---|---|---|
| Porte 1 | Non | Oui | Oui |
| Porte 2 | Non | Non | Oui |
| Porte 3 | Non | Oui | Non |
from theano import tensor as tt
with monty_model:
# Probability that Monty open each door
p_open = pm.Deterministic('p_open',
tt.switch(tt.eq(prize, 0),
# it is behind the first door
np.array([0., 0.5, 0.5]),
tt.switch(tt.eq(prize, 1),
# it is behind the second door
np.array([0., 0., 1.]),
# it is behind the third door
np.array([0., 1., 0.]))))
Monty ouvre la porte 3 contenant une chèvre
with monty_model:
opened = pm.Categorical('opened', p_open, observed=2)
Doit-on changer de porte?
with monty_model:
monty_trace = pm.sample(draws=10000, random_seed=SEED)
monty_df = pm.trace_to_dataframe(monty_trace)
fig
Donc, oui on devrait changer de porte
En option
We can also resolve the Monty Hall problem with pen and paper, as follows.
Throughout this calculation, all probabilities assume that we have initially chosen the first door. By Bayes' Theorem, the probability that the sportscar is behind door one given that Monty opened door three is
$$P(\textrm{Behind door one}\ |\ \textrm{Opened door three}) = \frac{P(\textrm{Opened door three}\ |\ \textrm{Behind door one})P(\textrm{Behind door one})}{P(\textrm{Opened door three})}.$$
The a priori probability that the prize is behind any of the doors is one third. From the table above, $P(\textrm{Opened door three}\ |\ \textrm{Behind door one}) = \frac{1}{2}$. We calculate $P(\textrm{Opened door three})$ using the law of total probability as follows:
$$ \begin{align*} P(\textrm{Opened door three}) & = P(\textrm{Opened door three}\ |\ \textrm{Behind door one})P(\textrm{Behind door one}) \\ & + P(\textrm{Opened door three}\ |\ \textrm{Behind door two})P(\textrm{Behind door two}) \\ & + P(\textrm{Opened door three}\ |\ \textrm{Behind door three})P(\textrm{Behind door three}) \\ & = \frac{1}{2} \cdot \frac{1}{3} + 1 \cdot \frac{1}{3} + 0 \cdot \frac{1}{3} \\ & = \frac{1}{2}. \end{align*} $$
Therefore
$$ \begin{align*} P(\textrm{Behind door one}\ |\ \textrm{Opened door three}) & = \frac{P(\textrm{Opened door three}\ |\ \textrm{Behind door one})P(\textrm{Behind door one})}{P(\textrm{Opened door three})} \\ & = \frac{\frac{1}{2} \cdot \frac{1}{3}}{\frac{1}{3}} \\ & = \frac{1}{3}. \end{align*}$$
Since $P(\textrm{Behind door three}\ |\ \textrm{Opened door three}) = 0$, because Monty wants the contestant's choice to be suspensful, $P(\textrm{Behind door two}\ |\ \textrm{Opened door three}) = \frac{2}{3}$. Therefore it is correct to switch doors, confirming our computational results.
Étude de cas: Sleep Deprivation¶
Les données proviennent de la librairie lme4 du language R, qui cite
Gregory Belenky, Nancy J. Wesensten, David R. Thorne, Maria L. Thomas, Helen C. Sing, Daniel P. Redmond, Michael B. Russo and Thomas J. Balkin (2003) Patterns of performance degradation and restoration during sleep restriction and subsequent recovery: a sleep dose-response study. Journal of Sleep Research 12, 1–12.
sleep_df.head()
Dans cette étude, chacun des sujets obtient son temps normal de sommeil lors du premier jour. Leur temps de sommeil est réduit à 3 heures pour les journées subséquentes. La colonnereaction est leur temps de réponse moyen à un nombre de tests. La colonne reaction_std est le résultat de la standardisation (Z-score) reaction sur tous les sujets et tous les jours.
fig
Chaque sujet devrait avoir un temps de réaction de base qui ne devrait pas trop différer les uns des autres
with pm.Model() as sleep_model:
μ_α = pm.Normal('μ_α', 0., 5.)
σ_α = pm.HalfCauchy('σ_α', 5.)
α = pm.Normal('α', μ_α, σ_α, shape=n_subjects)
Le taux de croissance des temps de réponse en fonction du nombre de jours de privation de sommeil différe selon chacun des sujets.
with sleep_model:
μ_β = pm.Normal('μ_β', 0., 5.)
σ_β = pm.HalfCauchy('σ_β', 5.)
β = pm.Normal('β', μ_β, σ_β, shape=n_subjects)
La combinaison du temps de réponse de base et du taux de croissance correspondent à nos observations.
with sleep_model:
μ = pm.Deterministic('μ', α[subject_ix] + β[subject_ix] * days)
σ = pm.HalfCauchy('σ', 5.)
obs = pm.Normal('obs', μ, σ, observed=reaction_std)
Ce type de model est connu sous le nom de model linéaire hiérarchique (mixe) "hierarchical (or mixed) linear model", parce que il permet à la pente et à la valeur initiale de variée selon le sujet, mais ajoute un terme de régularisation grâce à un apriori partagé par les sujets. Pour plus d'informations lire Gelman and Hill Data Analysis Using Regression and Multilevel/Hierarchical Models.
N_JOBS = 3
JOB_SEEDS = [SEED + i for i in range(N_JOBS)]
with sleep_model:
sleep_trace = pm.sample(njobs=N_JOBS, random_seed=JOB_SEEDS)
Diagnostic de la convergence¶
PyMC3 met à disposition une suite d'outils de diagnostic de convergence statistcal convergence diagnostics pour s'assurer que les échantillions soient une bonne approximation de la vrai distribution postérieure.
Les tracés d'énergie et la fraction bayésienne d'informations manquantes (BFMI) sont deux diagnostics pertinents pour l'échantillonnage Monte Carlo Hamiltonien, que nous avons utilisé dans cet exemple. Si les deux distributions dans le tracé d'énergie diffèrent de façon significative (surtout dans les queues), l'échantillonnage n'a pas été très efficace. La fraction bayésienne des informations manquantes quantifie cette différence avec un nombre compris entre zéro et un. Un BFMI proche de 1 est préférable, et un BFMI inférieur à 0,2 indique des problèmes d'efficacité.
Pour plus d'information : Michael Betancourt A Conceptual Introduction to Hamiltonian Monte Carlo et Robust Statistical Workflow with PyStan.
ax = pm.energyplot(sleep_trace, legend=False)
ax.set_title("BFMI = {:.3f}".format(pm.bfmi(sleep_trace)));
Contrairement aux tracés énergétiques et aux BFMI, qui sont spécifiques aux algorithmes de Monte Carlo Hamiltonien, la statistique Gelman-Rubin est applicable à n'importe quel algorithme MCMC, à condition que plusieurs chaînes aient été échantillonnées. Une statistique de Gelman-Rubin près de 1 est préférable, et les valeurs inférieures à 1,1 sont généralement considérées comme une indication de convergence.
max(np.max(gr_values) for gr_values in pm.gelman_rubin(sleep_trace).values())
Prediction¶
En inférence Bayesienne les prédictions sont effectuées en échantillionnant la distribution postiérieur, qui est la distribution des observations futures possibles sachant les observations antérieures.
with sleep_model:
pp_sleep_trace = pm.sample_ppc(sleep_trace)
pp_df.head()
grid.fig
Étude de cas: EEVEE¶
Statitique de l'échantillon¶
df.groupby('participant').mean().age_participant.plot.hist()
plt.title('age des paticipants')
plt.show()
df.astype(float).groupby('participant').mean().describe().iloc[:, 4:]
Statistique des avatars¶
show_avatar()
Concentrons-nous sur le réalisme¶
df_realisme = df.dropna(subset=['realisme'])
df_realisme[['realisme', 'avatar']].groupby('avatar').describe()
Premier model simple¶
On crée notre model génératif
Supposons que chaque observation du réalisme pour chacun des avatars provient d'une distribution normale, ayant une moyenne et une variance.
with pm.Model() as model_simple_realisme:
# realisme
## priors
μ_avatars_realisme = pm.Uniform('μ_avatars_realisme', 0, 100, shape=num_avatar)
σ_avatars_realisme = pm.HalfCauchy('σ_avatars_realisme', 25, shape=num_avatar)
## obsevations
obs_realisme = pm.Normal('obs_realisme',
mu=μ_avatars_realisme[df_realisme.avatar.cat.codes],
sd=σ_avatars_realisme[df_realisme.avatar.cat.codes],
observed=df_realisme.realisme)
display(model_to_graphviz(model_simple_realisme))
with model_simple_realisme:
trace_simple_real = pm.sample(draws=2000, cores=cores, chains=1)
pm.plots.traceplot(trace_simple_real)
plt.show()
pm.plot_posterior(trace_simple_real, varnames=['μ_avatars_realisme'])
plt.show()
pm.plots.forestplot(trace_simple_real, varnames=['μ_avatars_realisme'])
plt.show()
pm.plots.forestplot(trace_simple_real, varnames=['σ_avatars_realisme'])
plt.show()
On peut maintenant critiquer notre modèle¶
fig
On remarque que nos distributions ne concordent pas tellement lorsque les observations sont près de 100¶
On sous-estime donc le réalisme des avatars les plus réalistes
Utilisons une distribution Normale tronquée¶
entre 0 et 100
with pm.Model() as model_simple_realisme_tronquee:
# realisme
## priors
μ_avatars_realisme = pm.Uniform('μ_avatars_realisme', 0, 100, shape=num_avatar)
σ_avatars_realisme = pm.HalfNormal('σ_avatars_realisme', 25, shape=num_avatar)
## obsevations
obs_realisme = pm.TruncatedNormal('obs_realisme',
mu=μ_avatars_realisme[df_realisme.avatar.cat.codes],
sd=σ_avatars_realisme[df_realisme.avatar.cat.codes],
upper=100,
lower=0,
observed=df_realisme.realisme)
display(model_to_graphviz(model_simple_realisme_tronquee))
pm.plots.traceplot(trace_model_simple_realisme_tronquee)
plt.show()
On a permis aux moyennes d'être près des limites
pm.plots.forestplot(trace_model_simple_realisme_tronquee, varnames=['μ_avatars_realisme'], chain_spacing=0)
plt.show()
On peut maintenant comparer avec l'ancien modèle¶
fig
Comparaison avec le Widely Available Information Criterion
model_simple_realisme_tronquee.name = 'Normale Tronquée'
model_simple_realisme.name = 'Normale'
pm.compare({model_simple_realisme_tronquee:trace_model_simple_realisme_tronquee,
model_simple_realisme:trace_simple_real})
On peut à présent analyser nos inférence pour répondre à notre question:¶
Quels sont les avatars les plus réalistes?
pd.DataFrame(data={'realisme':trace_model_simple_realisme_tronquee['μ_avatars_realisme'].mean(axis=0), 'avatar':np.arange(1,17)}).sort_values('realisme', ascending=False).set_index('avatar').head()
Par contre, ces inférences sont affectés par le biais de notre échantillon, e.g. majoritairement des femmes dans la vingtaine.
Aussi, certains participants peuvent avoir un fort biais ayant un grand impact sur nos réponses.
Essayons de régler ces problèmes
Bâtissons un modèle hierarchique¶
Permettons à chacun des participants d'être biaisés.
On veut aussi estimer quelle est la variance dans le biais du groupe de participants
display(model_to_graphviz(model_herarchical_real))
pm.plots.forestplot(trace_model_herarchical_real, varnames=['μ_avatars_realisme'], rhat=False, chain_spacing=0)
plt.show()
pm.plots.forestplot(trace_model_herarchical_real, varnames=['μ_participant_bias_realisme'], rhat=False, chain_spacing=0)
plt.yticks(plt.yticks()[0], ['' if idx%4 else y for (idx, y) in enumerate(plt.yticks()[1]) ])
plt.show()
On peut toujours comparer avec les autres modèles
pm.densityplot([trace_model_simple_realisme_tronquee, trace_model_herarchical_real], models=['simple', 'herarchique'])
plt.show()
model_simple_realisme_tronquee.name = 'Normale Tronquée'
model_simple_realisme.name = 'Normale'
model_herarchical_real.name = 'Herarchique'
pm.compare({model_simple_realisme_tronquee:trace_model_simple_realisme_tronquee,
model_simple_realisme:trace_simple_real,
model_herarchical_real:trace_model_herarchical_real})
Modélisons maintenant en incluant plusieurs prédicteurs¶
display(graph)
Les prédicteurs liés au participant¶
pm.plot_posterior(trace_model_herarchical_real_all_bias_par, varnames=['β_age_par', 'β_ismale_participant', 'β_videogame_par', 'β_anime_par', 'β_travail_par',])
plt.show()
Les prédicteurs liés à l'avatar¶
pm.plot_posterior(trace_model_herarchical_real_all_bias_par, varnames=['γ_avatar_par_same_genre', 'β_avatar_genre', 'β_avatar_age',])
plt.show()
pm.plots.forestplot(trace_model_herarchical_real_all_bias_par, varnames=['μ_avatars_realisme'], rhat=False, chain_spacing=0)
plt.show()
On peut toujours comparer avec les autres modèles
model_simple_realisme_tronquee.name = 'Normale Tronquée'
model_simple_realisme.name = 'Normale'
model_herarchical_real.name = 'Herarchique'
model_herarchical_real_all_bias_par.name = 'Herarchique + covariables'
pm.compare({model_simple_realisme_tronquee:trace_model_simple_realisme_tronquee,
model_simple_realisme:trace_simple_real,
model_herarchical_real:trace_model_herarchical_real,
model_herarchical_real_all_bias_par:trace_model_herarchical_real_all_bias_par,
})
Simulation¶
En plus d'interpréter les données et de mieux comprendre les interactions entre le modèle et les différentes vairable, on peut utiliser notre modèle pour générer des données et simuler un nouvel échantillion, tout en gardant la notion d'incertitude/plausibilité.
Simulons un échantillon provennant d'une population plus diverse que celle obtenue par les liste de courriels de l'Université¶
fig
Finalement on peut comparer quels sont les avatars les plus réalistes en fonction du modèle¶
pd.concat([df_old.reset_index(), df_new.reset_index()] ,axis=1).head(8).style.hide_index()
Pour plus d'info¶
Statistiques bayésiennes¶
|
|
|
|
Plusieurs sont disponible à la bibliothèque ou en ligne avec Ariane 2.0
Échosystème de programmation probabiliste¶
| Probabilistic Programming System | Language | License | Discrete Variable Support | Automatic Differentiation/Hamiltonian Monte Carlo | Variational Inference |
|---|---|---|---|---|---|
| PyMC3 | Python | Apache V2 | Yes | Yes | Yes |
| Stan | C++, R, Python, ... | BSD 3-clause | No | Yes | Yes |
| Edward | Python, ... | Apache V2 | Yes | Yes | Yes |
| BUGS | Standalone program, R | GPL V2 | Yes | No | No |
| JAGS | Standalone program, R | GPL V2 | Yes | No | No |
| SPSS | Standalone program | Proprietary | ? | ? | ? |
Pas seulement pour les ingénieurs ou les statisticiens!¶
Plusieurs des développeurs de PyMC3 et Stan proviennent des sciences sociales, de psychologie et neurosciences.
Merci!¶
Le Jupyter notebook utilisé pour ces transparents est disponible https://github.com/alexisfcote/presentation_pymc3.